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1.  INTRODUCTION 


The  prime  practical  limitation  encountered  in  attempts  to  restore 
degraded  imagery  arises  from  noise  inherent  in  the  detected  image  data. 

The  most  basic  source  of  such  noise  lies  in  the  photon  fluctuations 
associated  with  the  detection  of  the  finite  amount  of  light  energy 
available  to  the  imaging  system.  Thus  photon  fluctuations  pose  a funda- 
mental limitation  to  the  "restorabi 1 i ty"  of  a degraded  image. 

In  the  first  part  of  this  report  we  develop  a model  which  can  be  used 
to  mathematically  and  statistically  describe  an  image  detected  at  low 
light  levels.  This  model  serves  to  clarify  some  of  the  basic  properties 
of  photon  noise,  and  provides  a basis  for  the  analysis  of  image  restoration 

to  fol low. 

In  the  second  part  of  the  report  we  consider  the  problem  of  linear 
least-square  restoration  of  imagery  limited  by  photon  noise.  The  form 
of  the  Invariant  least-square  restoration  filter  is  derived  using  the 
statistical  model  appropriate  for  photon  noise.  The  mean-square  error 
achievable  with  such  a filter  is  then  derived  in  terms  of  the  total 
number  of  photo-events  detected  in  the  image,  the  complexity  of  the 
object,  and  the  type  and  severity  of  the  image  blur.  Finally,  in  the  last 
part  of  the  report,  examples  are  presented  for  the  case  of  atmospherically 

degraded  images. 

2 . MnnPI IMG  OF  PHOTON-LIMITED  IMAGERY 

Our  model  for  the  detected  imagery  is  essentially  a two-dimensional 
analog  of  the  semi -classical  model  developed  by  Mandel  [1,2]  for  the  study 
of  photon  counting  statistics.  This  semi -class! cal  model  is  known  to 

yield  results  that  are  in  complete  agreement  with  a rigorous  quantum 


1 


mechanical  model  for  all  detection  problems  involving  the  photoelectric 
effect  [3].  Thus  the  vast  majority  of  image  detection  prob’ems  are 
properly  included  in  this  framework. 

2. 1 SEMI -CLASS  I CAL  MODEL 

The  model  utilized  is  that  of  an  inhomogeneous  or  compound  Poisson 
impulse  process.  Thus  the  detected  data  d(x,y)  is  represented  by 

N 

d(x,y)  = I 5(x-x^,y-y^)  (1) 

n=l 

where  6(-,-)  is  a two-dimensional  Dirac  delta  function, 
represents  the  location  of  the  n^^  photoevent  (i.e.,  the  release  of  a 
photoelectron),  and  N is  the  total  number  of  photoevents  produced  by 
the  image.  In  this  representation,  N , x^  and  are  all  regarded 

as  random  variables,  with  statistical  properties  to  be  described  in  the 
following. 

In  accord  with  the  semi -c  lass  i ca  1 theory  of  photodetection,  th.,- 
probability  that  N events  occur  in  an  area  A on  the  detector  is  taken 
to  be  Poisson, 


• 

A (x,y)dxdy 
J 

N 

P^(N)  = - 

.A 

— 0 V r\ 

X (x , y ) dxdy 
J J 

La  -j 

( 

n; 

exp 

where  the  "rate"  X(x,y)  is  related  to  the  classical  image  intensity 
i (x,y)  falling  on  the  detector  through 

X(x,y)  = , (3) 

hv 
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Here  n is  the  quantum  efficiency  of  the  photosurface  (assumed  independ- 


ent of  (x,y))  , h is  Planck's  constant,  v is  thp,  noan  optical 
frequency,  and  x is  the  detector  integration  time. 

Since  the  distribution  i (x,y)  of  classical  intensity  is  unknown 
a priori,  A(x,y)  must  ultimately  be  modeled  as  a random  process. 

However,  in  calculating  average  properties  of  the  detected  image,  it  is 
often  helpful  to  first  treat  A(x,y)  as  a given  known  function,  and  then 
average  the  results  over  the  statistics  of  A.  This  procedure  is  entirely 
consistent  with  Baye's  rule  of  statistics.  We  further  note  for  future  use 
that,  for  a given  A(x,y)  , the  "event  locations"  independent 

random  variables  for  different  n's  , with  common  probability  density 
function  [^] 


^(x  ,y  ) 
n n 

O)  ~ 

' > 

A (x,y)dxdy 

QO 


(M 


Note  further  that,  because  A(x,y)  is  proportional  to  the  classical 
intensity  , A(x,y)  ^0. 

Before  turning  to  an  examination  of  the  image  properties  implied  by 
this  model,  we  point  out  some  of  the  ways  it  can  be  generalized.  First, 
in  practice  the  photoevents  registered  by  a real  distributed  detector 
consist  of  finite  spatial  pulses,  not  delta  functions.  This  fact  can  be 
incorporated  in  our  model  by  passing  our  ideal  detected  data  (Eq.(l)) 
through  a linear  spatial  filter,  thus  spreading  the  delt-  functions  into 
pulses.  Second,  due  to  photoelectron  multiplication  noise,  the  areas  and 
shapes  of  the  various  pulses  may  themselves  be  random  variables.  This 
property  can  be  included  in  our  model  by  making  the  soatial  filter 
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randomly  space-variant.  These  sophistications  can  all  be  included  in  the 
model,  but  since  our  interest  lies  in  fundamental  limits,  there  is  little 
to  be  gained  by  incorporating  these  additional  non-fundamental  degradation 
We  close  this  section  by  showing  in  Fig.  1 a typical  classical 
intensity  distribution  and  a corresponding  typical  detected  image,  the 
illustration  being  one-dimensional  for  simplicity. 

2.2  SPECTRAL  DENSITY  OF  PHOTON-LIMITED  IMAGERY 

One  of  the  fundamental  properties  of  a detected  image  is  its  spectral 
density.  Our  goal  in  this  section  is  to  calculate  the  spectra)  density 
of  the  detected  image  described  by  Eq.(l).  If  D(vy,v  ) represents  the 
Fourier  transform  of  d(x,y)  , i.e., 


D(Vx.Vy) 

then  our  goal  is  to  calculate 


-]2t.  {'j  x+\>  y) 
d(x,y)e  dxdy 


(5) 


♦^,(vx,Vy)  = E[|D(vx,Vy)!^]  , (6) 

where  the  symbol  E[-]  indicates  an  expectation  over  the  statistics  of 

N , (x„,y^)  , and  X.  Before  makina  this  calculation,  we  must  first 

n n 

mention  a physical  restriction  we  impose  on  the  random  process  X(x,y)  , 

i.e.,  on  the  statistical  properties  of  the  classical  intensity  in  the 

image  plane.  * 

Any  image  produced  by  an  optical  system  must  contain  finite  total 

energy.  From  this  fact  it  follows  that  for  every  sample  function  of 

the  random  process  X(x,y)  , 

00 

[[  X{x,y)dxdy  < » . (7) 
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This  equation  is  sufficient  to  imply  that  any  particular  sample  function 
^(x,y)  has  a well-defined  Fourier  transform. 


A(Vx»^y) 


J 


X (x,y)e 


(8) 


Our  goal  now  is  to  relate  the  spectral  density  ‘*’(j  detected 

image  to  the  spectral  density 


,2 


(9) 


where  the  expectation  is  over  the  statistics  of  \(x,y). 

To  calculate  the  spectral  density  detected  image, 

we  first  Fourier  transform  Eq.(l),  with  the  result 


N ■J2'(v  X + y ) 

#<fc/  \ K'  AnYn 

0(^'x*'-’y)  = I e 


(10) 


n=  I 


The  squared  modulus  of  this  quantity  is 


|0(vx,Vy) 1 


N N 

I I e 

n= I m= 1 


■JZtt 


X n m Y n m 


. (11) 


It  remains  to  find  the  expected  value  of  this  quantity  over  the  statistics 

of  N , X(x,y).  It  is  convenient  to  first  yegard  N and 

X(x,y)  €s  given  (known)  quantities,  average  over  the  statistics  of 

(x  ,y  ) and  (x  ,y  ) , and  then  average  over  N and  X.  Thus  our  first 
' n ^n  m’ ^m  ^ 

goal  is  to  compute 


E 

nm 


|D(v^. 

N 

I 

n=l 


v„(x  -X  )+v„(y  -y  ) 
X n m Y'^n  'm 


(12) 


- 5 - 


Wlfn 


where  E signifies  an  average  over  (x  ty  ) and  (x  ,y  ). 
ntn  n n m rn 


Tvx>  classes  of  terms  can  be  identified.  First,  there  are  N terms 
for  which  n * m , each  of  which  yields  unity.  Second,  there  are 


*■  N terms  for  which  n d m.  For  such  terms  we  kno^j  that  (x  ,y  ) 

n n 


and  (x  ,y  ) are  independent  random  variables, and  therefore  that 
ro  m 


A(x  ,y  ) 

n n m m 


ty-) 
m m 


03) 


jj  A(x,y)dxdy  X(x,y)dxdy 


For  these  * N terms,  the  result  of  the  averaging  process  is 


run 


“j2n (v^x+Vyy) 
X(x,y)e  dxdy 


X(x,y)dxdy 


(H) 


The  result  of  averaging  |0(vj^,Vy)l  over  the  statistics  of  (^„iy„) 


n 'n 


and  (x  ,y  ) becomes 
ro  m 


A ( Vx  • ^ 


A(0,0) 


(15) 


Continuing  our  averaging  process,  we  next  find  the  expected  value  of 
over  the  random  variable  N , given  X(x,y).  Representing  the 


conditional  mean  of  N (given  X)  by  , and  noting  that,  for 


Poisson  statistics, 
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:["'-•<]  = (R(j)) 


= jj  A(x,y)dxdy  = A(0,0) 


we  find  that 


'n,m,N  1^]  " '^(A)  - (17) 


Finally,  averaging  over  the  statistics  of  A.(x,y)  , we  obtai 


>>d(vx.VY)  = N + 4>^(vj^.Vy) 


where  N is  the  unconditional  mean  of  N. 


Thus  the  spectral  density  of  the  detected  image  is  the  sum  of  a 


4)^.  Alternate  forms  of  this  result  are  also  useful.  First,  if  we  define 


a normalized  spectral  density 


^x(Vx.Vy) 


^/Vx.Vy) 


♦^(o.o) 


then  we  have  that 


*d^''x*''y^  * N + (v^. Vy) 


Furthermore,  since  A(x,y)  is  proportional  to  the  classical  intensity 
i (x,y)  , we  must  have 


where  ♦j(vjj.Vy)  Is  the  spectral  density  of  the  classical  intensity, 
normalized  to  unity  at  '^x  **  '^Y  * Thus  we  write  our  final  result, 

*d^'"X'''Y^  " N + (N)^4j(vx,vy)  . a 


constant  spectral  level  N , plus  the  spectral  density  cf  the  rate  function, 
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We  illustrate  this  result  in  Fig.  2. 

2.3  FLUCTUATIONS  OF  SPECTRAL  DENSITY  ESTIMATES 


i 


In  a certain  class  of  imaging  problems,  the  desired  end  result  is  an 
accurate  estimate  of  the  normalized  spectral  density  <^  . of  the  classical 
image  intensity.  For  example,  such  is  the  case  for  Labeyrie's  speckle 
interferometer  [5].  which  is  currently  of  great  interest  in  astronomy. 
Because  of  the  simple  relation  that  exists  between  <I>.  and  the  spectral 
density  of  the  detected  image,  a reasonable  approach  is  to  first 

estimate  4>^  , then  express  <t>.  as  (c.f.  £q.(22)) 


(N)^ 


(23) 


The  quantity  N is  simply  a measure  of  the  total  image  brightness,  which 
we  assume  Is  either  known  a priori  or  can  be  determined  accurately  by  a 
suitable  photometric  measurement.  Thus  the  fluctuations  in  our  estimate 

A 

of  <j>j  will  be  determined  by  the  fluctuations  inherent  in  our  measurement 

of  It  Is  these  fluctuations  that  we  wish  to  find  here. 

2 

An  estimate  of  is  made  by  measuring  |0(vjj,Vy)|  for  a single 

detected  image.  The  expected  value  of  this  measurement  is,  of  course, 
<lid(vx,Vy).  But  how  far  off  from  the  expected  value  is  a single  measure- 
ment likely  tc  be?  To  answer  this  question,  we  must  find  the  second 
moment  of  |0|  , l.e.,  we  must  calculate 

r N N N N r ( r 

e[|d1]  » I I I I ^[exp  -J2Tr[vx(x^-x^  H-  X -X  ) 

n»l  m«l  p"l  q“l  U ' 

vVjl]  • 

This  calculation  is  a lengthy  one  and  is  deferred  to  the  Appendix. 

For  a known  (non-random)  A(x,y)  , the  result  Is 
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It  iww  II  liill■lll^> 


e[|D|'‘]  - N(^)  + 2(n,^,)2  + n(l+N,j))|A(v,.v^)|^ 

+ |a(2vjj,2vy) 

. 4-  A(2Vj,,2vY)[A*(Ojj,VY)j^  + A*(2Vj^,2Vy)  j^A(vjj,Vy)J 

k 

+ lA(Vj^,Vy)  I . 

Representing  A(vy,Vy)  in  terms  of  a modulus  and  phase. 


A(Vy,Vy)  = jA( 


Vy.Vy) |el 


0 ♦ ^y) 


we  can  equivalently  express  the  second  moment  of  |d(  as 

e[|o|'']  ■ V)  *e(V))^*‘'('^"(a))I'“''x’VI^ 

+ |A(2\)jj,2vy) 

+ 2|A(2Vj^,2vy)  I |a(vj^,Vy)  l^cos[e(2Vj^,2Vy)  - 2e(vj^,VY)] 

+ 1A(Vj(.Vy)|^  . (27) 

The  task  now  remains  of  averaging  over  the  statistics  of  X(x,y). 

If  the  image  intensity  distribution  extends  over  a region  of  size  LxL  , 

then  under  rather  general  conditions,  the  central  limit  theorem  can  be  used 

to  show  that,  for  ^ and  Vy  » I"  . A(vj(,Vy)  Is  approximately  a 

circular  complex  gaussian  random  process,  with  correlation  extending 

over  a region  of  dimensions  ^ l’  frequency  domain.  It  follows 

2 

that  the  phase  6 is  uniformly  distributed  on  (-it,TT)  , and  that  |a| 
obeys  negative  exponential  statistics.  Furthermore,  for  such  frequencies 
e(2vjj,2vj()  , e(Vy,Vy)  , lA(2vy,2vy)i  and  |A(vjj,Vy)|  are  all  approximately 


.JJ,-Vy 


independent.  From  these  facts  we  conclude  that 


.•ij>«.T;i\flrT.^;^W^'?iW<T'’rrf^ 

" P»ttW>^'‘V'‘'*''n‘i'«^'^ 


e[|a(vjj,Vy)  l^j  “ 

E j^|A(2vjj,2vy)  l^j  ■ 4>^  (2vjj,2vy) 

e[21a(2vjj»2\>y)  1 |A(Vjj,Vy)  j^cosje(2\)j^,2vY)  “ 2e(vjj,Vy)]j  = 0 

(28) 

2 


E[lA(Vjj,Vy)  1^]  - 24>^(Vj^,Vy) 


Now  If  we  subtract  the  square  of  the  mean  of  |0 
of  Eq.(l8),  we  obtain  the  variance  of  |d|  , 

2 * ‘ N + (N)^  + 2(2+N)$^(vj^,Vy) 


I.e.,  the  square 


+ ♦^(2Vj^.2Vy)  + »;,(Vj(,Vy) 


(29) 


Equivalently,  using  (21),  we  have 

i2 


N + (N)^  + 2(24.N)(N)^*j(vj^,Vy) 


+ (N)^*?!  (2vjj,2vy)  + (N)^4>^(vjj,Vy) 


(30/ 


An  Important  conclusion  can  be  drawn  directly  from  Eq.(30).  The 
fluctuations  of  the  spectral  density  of  the  detected  Image  at  frequency 
(vj^,Vy)  depend  not  only  on  the  spectral  density  of  the.  classical  intensity 
at  (v^,Vy)  , but  also  on  that  spectral  density  at  frequency  (2vj^,2vy)  ! 
Stated  In  other  words,  a frequency  component  of  the  classical  Intensity 
at  frequency  (2v^,2vy)  Induces  fluctuations  of  the  spectrum  at 
(vx»Vy).  This  "ha If- frequency  noise"  phenomenon  is  a fundamental  property 
of  photon  limited  Images.  It  has  been  noted  previously  in  a different  but 
related  context  by  Walkup  [6]. 

We  close  this  section  by  presenting  an  expression  for  the  rms 
signal-to-noise  ratio  associated  with  an  estimate  of  at  frequency 
(vxtVy).  Subtracting  the  bias  N associated  with  the  mean  of  |0|  , 
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we  obtain 


eLIdI^I  - N 


J (Vx.Vy) 


}i?(Vx.Vy)  + 


i 


Note  that  as  N > “ , the  r.m.s.  signal-to-noi se  ratio  approaches  unity, 
in  agreement  with  the  classical  results  on  the  statistical  fluctuations 
of  per lodograms  [7]  • 

Figure  3 illustrates  the  conclusions  of  this  analysis  with  a specific 

•example.  The  classical  intensity  distribution  is  taken  to  be  a sinusoidal 

fringe  of  frequency  Vq  and  length  L.  A typical  sample  function  of 
2 

|D(v)1  is  shown.  Note  that  excess  fluctuations  are  present  at  Vq/2 
due  to  the  presence  of  the  fringe  at  frequency  Vq. 

3 - LINEAR  LEAST-SQUARES  RESTORATION  OF  DEGRADED 

PHOTON-LIMITED  IMAGES 

In  many  practical  problems  of  interest,  the  detected  image  data  arises 
from  a blurred  image  of  the  object  of  interest.  For  example,  the  object 
may  suffer  significant  motion  during  the  detection  interval  i , thus 
blurring  the  image.  Alternatively,  and  of  greater  interest  here,  the 
detected  image  may  be  seriously  degraded  by  the  spatial  and  temporal 
fluctuations  of  the  refractive  index  of  the  atmosphere,  i.e.,  by  "atmospheric 
seeing".  At  low  light  levels,  the  detected  Image  further  suffers  from 
photon  noise  of  the  type  discussed  in  the  previous  section.  In  order  to 
extract  as  much  information  as  possible  about  the  object  from  the  detected 
image  data,  we  seek  a method  of  image  restoration  which  will  enhance  object 


II 


detail  without  unduly  emphasizing  the  noise  assoc iated  with  the  discrete 
photo->events  composing  the  detected  Image.  In  the  sections  to  follow, 
we  consider  one  approach  to  this  problem. 

3.1  LINEAR  LEAST-SQUARES  RESTORATION 

The  approach  we  shall  Investigate  here  is  commonly  referred  to  as 
I inear  least-squares  restoration.  The  philosophy  behind  this  approach  is 
perhaps  best  explained  with  the  aid  of  Fig.  k.  The  function  o(x,y) 
represents  the  true  object  brightness  distribution,  or  alternatively  the 
image  that  would  be  produced  by  an  ideal  optica]  system  (free  from 
aberrations  and  free  from  any  blur  due  to  diffraction)  and  an  ideal 
noise-free  detector.  Our  goal  Is  to  determine  o(x,y)  from  the  actual 
detected  data  with  the  greatest  possible  accuracy.  Following  the  upper 
branch  of  Fig.  k,  the  ideal  object  suffers  a perfectly  known  blur  on  passage 
to  the  image  plane,  this  blur  being  Introduced  by  diffraction,  fixed 
aberrations,  and  other  external  causes,  such  as  object  motion,  atmospheric 
seeing,  etc.  We  assume  that  all  of  these  blurs  can  be  lumped  together 
and  represented  by  a single  known  linear,  space- invariant  filter,  with 

A 

impulse  response  b(x,y)  , or  optical  transfer  function  B(vj^,Vy), 


B(Vjj,Vy) 


II  b(x,y) 


-j2u(v  x+v  y) 
e * ’ dxdy 


•00 


II 

•00 


b(x,y)dxdy 


(32) 


To  represent  the  statistical  fluctuations  introduced  by  the  detection 
process,  the  blurred  image  is  now  applied  to  a "Poisson  generator",  which 
produces  a Poisson  impulse  process  with  rate  A(x,y)  proportional  to  the 
intensity  of  the  blurred  image.  The  output  of  the  Poisson  generator  Is 
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the  detected  Image  data  d(xty)  , on  which  we  must  base  our  estimate  of 
o(x,y). 

Our  restoration  procedure  is  to  apply  the  detected  image  data  to  a 
linear,  sp^ce- invar  rant  restoration  filter  with  impuise  response  h(x,y) 
and  optical  transfer  function  H(uj^,Vy).  This  transfer  function  will  be 
chosen  to  minimize  the  expected  value  of  the  mean-squared  error. 


_ « 


2 


e 


(x,y)dxdy 


L>oo 


{ r (x,y)-o(x,y) I dxdy 


(33,1 


where  the  error  e(x,y)  represents  the  difference  between  the  restored 
image  r(x,y)  and  a certain  filtered  version  of  the  object,  o(x,y)  , 
the  expectation  being  over  the  statistics  of  o(x,y)  and  the  statistics 
of  the  detection  process. 

The  choice  of  a filtered  object  o(x,y)  , rather  than  the  true  object 
o(x,y)  , for  defining  the  error  requires  some  comment.  The  restoration  of 
object  frequency  components  beyond  the  d i ffract ion- I imi ted  cutoff  of  the 
optical  system  is  impossible  to  achieve  with  any  linear  invariant 
restorat(on  filter.  Hence  the  best  we  can  hope  to  accomplish  is 
restoration  of  those  frequency  components  lying  within  the  diffraction- 
limited  passband.  Accordingly,  we  count  as  error  only  the  differences 
between  the  restored  spectrum  R(vjj,Vy)  and  the  (possibly  modified) 
portion  of  the  spectrum  lying  within  the  observable  passband.  The  frequency 
spectrum  of  o(xty)  is  thus 


t'(VjjtVy) 


Vy  S(Vjj,Vy) 
0 


(Vx.Vy) 

in  observable  passband 
otherwise 


- 13  - 


Usually  we  will  take  S(vjj»Vy)  to  be  of  the  form 


S(VjJ,Vy) 


(vj^,Vy)  In  observable 
passband 


otherwl se 


However,  we  note  that  the  resulting  o(x,y)  can  have  negative  values  in 
this  case.  We  could  alternatively  choose  S(vj^,Vy)  to  be  the  diffractic 
limited  optical  transfer  function  of  the  system,  thus  guaranteeing  a 
positive  o(x,y).  In  the  analysis  to  follow  we  leave  S(vj^t>Vy)  completely 
general,  but  eventually  we  choose  the  form  of  Eq.(35)  for  mathematical 
simplicity. 

Some  final  comments  are  in  order  regarding  the  non-opt imal i ty 
of  the  restoration  procedure  described  above.  First,  it  is  well  known 
that  linear  least-squares  restoration  is  not  optimal  *n  the  sense  of 
maximum  likelihood  or  maximum  a posteriori  probability  when  the  image 
statistics  are  Poisson.  Rather,  non-linear  filtering  is  required  for  true 
optimality  i7] . Secondly,  the  choice  of  a space- invar iant  linear  restor- 
ation filter  undoubtedly  reduces  performance  even  further;  it  seems  clear 
intuitively  that  a space-variant  filter  can  perform  better  than  a space- 
invariant  filter  In  the  presence  of  signal-dependent  noise.  However,  it 
should  be  pointed  out  that  both  optimal  non-linear  filtering  and  optimal 
space-variant  linear  filtering  are  in  general  computationally  less 
efficient  than  linear  space- invar iant  filtering.  For  this  reason,  there 
remains  a strong  Interest  In  knowing  the  limitations  of  linear  space- 
invariant  least-squares  restoration  for  photon- 1 Im I ted  imagery. 
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3.2  THE  FORK  OF  THE  RESTORATtON  FILTER  AND  THE  QUALITY 
OF  THE  RESTORED  IMAGE 

In  this  section  we  first  derive  the  form  of  the  linear  space- 

invariant  least-square  restoration  filter  for  photon- 1 imi ted  images. 

Our  goal  is  to  choose  a filter  transfer  function  H which  minimizes 

2- 


//  e^(x,y)dxdyj.  By  Parseval's  theorem,  it  is  equivalent  to  minimize 


L-oo 


where  6 is  the  Fourier  transform  of  e. 


Interchanging  orders  of  integration  and  expectation,  we  find  that  it 
suffices  to  minimize  at  each 


[i6(^x''V^*^]  “ e[|dh  - o;^ 


(36) 


where  0(v^,Vy)  is  the  Fourier  transform  of  o(x,y).  The  m i n imi zat ion  is 
straightforward  and  yields 


H(Vx*'''y)  “ “ 


‘do 


'd 


(37) 


where  is  the  spectral  density  of  the  detected  image  d(x,y)  , while 

is  the  cross-spectral  density  of  d(x,y)  and  o(x,y)  , 


■^do 


(38) 


Straightforward  cal cu lat ions,  using  the  Poisson  impulse  model,  show  that 
‘t'^(vx.Vy)  = N + (N)^| BCv^t Vy)  I ^'5'^ ( , vy ) 


‘^do^^X’  Y^  “ S{vx.Vy)B  X’^Y^ 


(39) 


We  conclude  that  the  transfer  function  of  the  restoration  filter  is  given 
by 
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H(Vjj,\Jy) 


N S(vj^,Vy)B  ('''.^iVY)'&Q{vj^,Vy) 

I •♦•  N|B  (vjj,Vy)  |^4>^(yj^,Vy) 


m 


Note  that  at  zero  spatial  frequency  the  gain  of  this  filter  is 
N/(I+N)..  Thus  the  normal  ized  restoration  transfer  function  is 


H(Vj^,Vy) 


(l+N)S(vjj,Vy)B  (^j{»yy)*Q(vjj,yy) 


1 + N|6(w^,yy)|  <t  (vj^,Vy) 

We  turn  next  to  calculation  of  the  total  mean>squared  error  c 
achieved  when  the  restoration  filter  of  Eq. (40)  is  used.  Using  (36) 
we  see  that 


(41) 


e ® E 


1 4(Vy,Vy)  I dvj^  dvy  I 
[|H|^^  - «<^do  ■ '^"‘^do  ^ 

•00 

Substituting  (39)  and  (4o)  in  (42),  we  obtain  after  some  algebra 


(42) 


"X  • 


(N)^|S|  V 


I * NIB 


(43) 


The  total  mean-square  error  In  the  Image  is  not  a particularly 
meaningful  quantity  in  itself.  Note  in  particular  that,  as  N grows 
large,  so  too  does  r , in  direct  proportion  to  N.  However,  the  mean- 
square  value  of  the  object  within  the  observable  passband  is 
(N)^  //  1^1  ^Qdvj^dvy  ; hence  the  "signal"  component  of  the  output  power 

“CO 

— 2 

rises  in  proportion  to  (N)  , yielding  a net  increase  in  restored  image 

quality  as  N increases. 
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It  is  highly  appealing  to  define  a single  parameter  Q to  represent 
the  overall  quality  of  the  restored  image.  Many  definitions  of  such  a 
parameter  are  possible,  but  none  can  be  fully  Justified  as  being  the 
best  conceivable  choice.  Here  we  shall  consider  two  possible  choices 
for  a definition  of  H. 

Our  first  definition  depends  on  the  mean-square  error  > , which  takes 
into  account  both  the  statistical  fluctuations  of  the  restored  image  due 
to  the  photodetection  process  and  the  defects  of  the  restored  image  caused 

O' 

by  residual  uncompensated  blur.  Noting  that  as  N 0 , t >•  (N)  //  [s] 
♦^dVjjdVy  , we  find  the  proper  definition  of  the  quality  factor  to  be 


(N)^ 


:,2: 


jj 


[S| 


or  equivalently 


Q 


I 


I |S|\dv^dVy 

-to 

•7  Isl^^pdv^dvy 

• I + iT1b| V 

00  * O 


A second  possible  quality  parameter  which  has  a great  deal  of  appeal 
will  now  be  discussed.  This  second  parameter  is  defined  to  be 

(I2  » (Av)^log(l  + (45) 

2 

where  (Av)  is  a measure  of  the  two-dimensional  restored  bandwidth, 
while  S/N  is  a measure  of  the  mean-square  signal-to-noise  ratio  in  the 
restored  image.  The  similarity  of  this  quality  measure  to  a measure  of 
information  content  is  apparent.  A reasonable  definition  of  the  restored 
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bandwidth  Is 


I*/ 

i 


t 


1: 


(Av)^  - II  |BH|dvjjdvy  (46) 

where 

(l+?a)S(V„,U„)  |B(v„,Vy)  1^*  (Vy.Vy) 

8H  ^ . (47) 

I ♦ N|B(Vjj,Vy)  I ♦o(''j^*Vy) 


As  for  the  parameter  S/N  , the  choice  of  a definition  is  less  clear 
cut.  One  suitable  choice  is  the  mean-square  signal-to-noise  ratio  per 
degree  of  freedom  of  the  restored  image. 


1 . N 

**  ’ A,(4v)^ 


(48) 


where  A|  represents  the  area  occupied  by  the  restored  image.  Only  the 
first  of  the  above  quality  measures  (Q|)  will  be  used  in  section  4. 

3.3  THE  OEPENOENCE  OF  THE  WORHALIZED  SPECTRAL  DENSITY 
OF  THE  OBJECT  ON  OBJECT  COMPLEXITY 

The  results  of  the  previous  sections  have  demonstrated  that  the 

degree  to  which  restoration  is  possible  depends  not  only  on  the  total 

light  flux  (H)  and  the  optical  transfer  function  of  the  blur  (B)  , but 

* 

also  on  the  normalized  spectral  density  of  the  object  (4>^)>  In  this 
section  we  explore  the  dependence  of  this  normalized  spectrum  on  "object 
complexity".  To  explore  this  question  in  a completely  general  way  is 
extremely  difficult.  Accordingly,  we  examine  two  specific  models  of  the 
object,  neither  of  which  is  entirely  realistic,  but  from  which  we  can 
deduce  trends  valid  for  more  general  object  models. 


I 
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A. 


For  ccMl^>ar{son  purposes  we  note  that  the  simplest  possible  object 
is  an  ideal  point  source,  with  brightness  distribution 

o(x,v)  - o^6(x-Xj,y-y,)  . (i»9) 

While  such  an  object  cannot  exist  physically,  nonetheless  it  serves  as 
a useful  idealized  case  against  which  we  can  compare  the  results  of. 
complicated  object  models.  For  the  point  source  described  above,  the 
normalized  spectral  density  is  given  by 

^’o^Vx-yy)  • I (all  ^x'^Y^  • 

The  first  model  utilized  to  represent  a more  complicated  object  is 
a natural  general izat ion  of  the  case  of  a single  point  source.  We 
suppose  that  there  exist  H equally  intense  point  sources, 

M 

o(x,y)  » I <>o'5{x-x^,y-y^)  . (51 ) 

ni»l 

We  suppose  that  the  locations  (x^,y^)  are  independent  random  variables 
uniformly  distributed  over  a square  object  field  of  size  L»L.  Omitting 
the  calculations,  which  are  straightforward,  we  find  the  normalized 
object  spoctral  density  in  this  case  to  be 

which  is  shown  in  Fig.  5^*  Note  that  an  increase  of  object  complexity 

has  led  to  a decrease  in  the  level  of  the  normalized  spectral  density, 

except  at  extremely  low  spatial  frequencies. 

What  is  the  effect  of  object  complexity,  represented  by  N , on  the 

quality  of  the  restored  image?  Using  the  fact  that  ■ — over 

OH 

most  spatial  frequencies  of  interest,  examination  of  Eqs.  (44),  (45) 

9 - 
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and  (48)  shows  that  Qj  and  primarily  functions  of  the  total 


number  of  photoevents  contributed  by  a single  point  source. 


_N^ 

M 


(53) 


If  n is  held  fixed,  Q is  essentially  independent  of  how  many  point 
sources  are  in  the  field. 

A second  and  somewhat  more  general  model  for  the  object  can  be  formu- 
lated as  follows.  Let  the  object  be  modeled  as  a stationary  randcwn 
process,  o(x,y)  0 , with  mean  o _ 0 and  variance  We  further 

suppose  that  this  stationary  randotn  process  is  confined  to  an  L L 
square,  and  is  zero  outside  that  square.  Ihis  space  limitation  can  be 
explicitly  introduced  by  multiplying  the  stationary  o(x,y)  by  a window 
function  rect (x/L) rect (y/L) . We  wish  to  calculate  the  normalized 
spectral  density  for  this  object  model.  The  calculation  is  tedious  but 
again  straightforward.  To  state  the  results  in  succinct  form,  we  define 

the  following  additional  symbols: 

2 

A • L represents  the  area  of  the  object  t ield; 

A^  represents  the  correlation  area  of  the 

random  process  o(x,y)  , and  is  specifically 
defined  as 


Jj 


»^(Ax,Ay)dAxdAy 


where  is  the  autocovar lance  of  o(x,y)  , 

normalized  to  unity  at  the  origin;  and 


♦ {v„,v,,)  is  the  spectral  density  of  the  fluctuations 

qO  a Y 

of  the  object  about  its  mean,  normalized  to  unity 
at  the  origin. 


20 
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With  these  definitions,  the  normalized  spectral  density  can  be  expressed 
approximately  as 


♦o(Vx.Vy) 


■m 

(■f)  X 


. 2 2 
sine  LVj^  sine  LVy 


(54) 


where  the  chief  approximation  is  that 


(55) 


Figure  5b  shows  a typical  plot  of  this  norma  1 ized  spect ra I density. 

Of  most  importance,  we  note  that,  except  at  the  lowest  spatial  frequencies, 
the  normalized  spectral  density  has  a value  less  than  (n^/o)^ (A^/A) . 

This  parameter  plays  a role  similar  to  (M)  ' in  the  previous  model.  In 
this  case  we  define  the  parameter  n as  the  mean  number  of  photoevents 
contributed  by  a single  correlation  area  of  the  object. 


Again  referring  to  Eqs.(44),  (45)  and  (48),  we  find  that  the  quality 

— - 2 

of  the  restored  image  will  depend  primarily  on  the  parameter  n(o^/o) 
for  any  given  blur. 

4.  APPLICATION  TO  ATHOSPHERtCALLY  DEGRADED  IMAGES 

Our  goal  here  is  to  apply  the  results  of  section  3 to  the  specific 
case  of  atmospherically  degraded  images.  All  imagery  considered  here 
will  be  assumed  to  be  recorded  with  an  exposure  time  that  is  much  longer 
than  the  characteristic  fluctuation  time  of  the  atmosphere.  Thus  we  arc 


dealing  only  with  long-exposure  imagery.  However,  we  consider  both 


ordinary  long-exposure  Imagery  and  "tilt-removed"  long-exposure  Imagery. 


In  the  second  case  It  is  assumed  that  a perfect  tilt  removal  system 


operates  to  keep  the  Image  perfectly  centered  on  a fixed  point  at  ail 


times. 


Using  the  results  of  Fried  [8],  we  have  that  the  OTF  of  atmospherically 


induced  blur  is  given  by 


B^(v)  » expj-3.A4(^)  1 - cc(Y-)  ^ i 


where 


is  the  mean  wavelength^ 


is  the  focal  length  of  the  telescope; 
represents  the  telescope  diameter; 
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/vw  + Vv  represents  radius  in  the  spatial  frequency 


plane; 


represents  the  coherence  parameter  of  the  atmospheric 


wavefront  distortions,  as  defined  by  Fried  [8];  and 


a takes  on  the  value  zero,  one  half,  or  one,  according  to 


whether  the  image  is  recorded  with  no  tilt  removal 


(«=0) ; with  tilt  removal  and  "far  field"  atmospheric 


propagation  conditions  [8]  i or  with  tilt  removal 


and  "near  field"  atmospheric  propagation  conditions 


[8]  (a=l). 


In  addition  to  the  atmospherically  Induced  blur,  we  assume  that  blur 


is  introduced  due  to  diffraction  by  the  finite  aperture  size  of  the 


telescope.  For  a perfect  circular  telescope,  we  have  an  optical  transfer 


function 


J 
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BjM 


2 -I 
~ cos 

TI 


(y-W'-(y 


for  V < V , zero  otherwise,  where  v = D/(AF)  represents  the 
— o o 

diffraction  limited  cutoff  frequency.  For  computational  purposes,  a 
convenient  approximation  due  to  Hufnagel  [9]  can  be  used, 


BtW  i I - 


The  total  OTF  of  the  imaging  system  is  simply  the  prcc'jct  of  Eqs.  (57) 
and  (58)  or  (59), 

B(v)  = B.p(v)B^(v)  . (6( 

Incorporating  the  definition  of  in  Eq.(57)»  we  obtain  for  the 


total  OTF, 


B(v);j^l  - l.25(-y  *0,25(^)"j 


for  V <_  , zero  otherwise. 

2 

We  have  used  numerical  integration  to  calculate  (Av)  (Eq.(46)), 
and  Qj  (Eq.(44))  for  the  case  of  a point-source  object  (l*^  =1)  with 
blur  6 of  Eq.(6l),  and  an  ideal  transfer  function  S of  Eq.(35). 

In  this  case. 


(Av)^  = 2n(l+N)  f 

L I + n|b(v) 


and  we  can  show 


(63) 


' TtV^  - (Av) 
o 

In  the  computations  we  have  assumed  the  parameter  values 


r = 10  cm  and  5 cm 
o 

D “ 152.5  cm 

_ _c 

X = 5x10  cm  . 

It  Is  not  necessary  to  assume  a specific  focal  length  for  the  telescope 

if  we  work  with  spatial  frequencies  Q ^ pv  measured  in  cycles  per 

radian  of  arc  In  the  sky.  The  ratio  of  in  Eq.(6l)  is  replaced 

2 2 2 

by  , where  , and  we  calculate  (Afi)  = F (Av)  rather 

2 

than  (Av)  . Likewise,  Q|  is  re-defined  as 


(An) 

- (An)^ 

O 


(64) 


with  no  change  in  its  numerical  values  resulting. 

in  Fig.  6 we  show  plots  of  the  "maximum  restorable  frequency" 

AnZ/rT  vs.  N for  the  case  of  a point-source  object  and  the  parameter 
values  specified  above.  The  incident  light  flux  is  varied  over  eight 
orders  of  magnitude.  Figure  6a  corresponds  to  the  "no  tilt  removal" 
case,  while  6b  and  6c  represent  the  "tilt  removed"  cases  for  far-field 
and  near-field  atmospheric  propagation  conditions,  respectively.  Note 
that  in  all  cases,  the  maximum  restorable  frequency  increases  very  slowly 
with  N , Implying  that,  for  the  parameter  values  specified  here,  truly 
enormous  amounts  of  light  flux  are  necessary  to  record  an  image  which 
can  be  substantially  restcred. 

To  illustrate  this  point  further,  consider  the  transfer  function 

A A 

BH  of  the  cascade  of  the  blur  '^nd  the  deblur.  From  Eq.(47)  with 
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N » J and  ♦ * 1 , we  see  that  wJthIn  the  diffract  ton-limited  passband 


we  have 


BH 


1 + n1b|^ 


(65) 


If  we  wish  to, restore  the  frequency  component  at  v to  Jg  the  amplitude 
It  had  before  the  blur,  we  require  that  N|B|  = 1 , or 

N = |8|"^  . (66) 

Assuming  that  a = 0 (no  tilt  reraovol)  and  for  simplicity  neglecting 


the  d If f ract Ion- 1 iml ted  portion  of  the  OTF,  we  require  from  £q.(6l)  that 

5/3.  v5/3^ 

(67) 


N = expj6.83(| 


I 


) • 


V,'ith  D “ 152  cm  , fQ  ® 10  cm  and  A = 5><10  ^ cm,  we  find 

if  - expje.as^Jj^  j 

with  expressed  In  cycles  per  mill  trad  ian.  Now  for  u = il^/2  = 1520 
cycles  per  mil  11  rad ian,  we  find  that 


•«  37 

N * 6x10  ^ photoevents  . 


.87 


Thus  more  than  10  ' photoevents  are  required  in  the  detection  process 
to  achieve  this  degree  of  restoration! 

In  Fig.  7 we  have  plotted  the  quality  factor  Qj  for  the  same 
range  of  N , again  for  two  values  of  r^.  Figures  7a,  b and  c 


liiiiiiiiaiiittli 
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again  correspond  to  no  tilt  removal,  tilt  removed  with  far-fleld  atmos- 
pheric propagation,  and  tilt  removed  with  near-fleld  atmospheric  propa- 
gation. Again,  for  the  parameters  chosen,  there  is  only,  a small  change 
of  image  quality  over  this  wide  range  of  N*. 

These  results  support  the  experimentally  observed  fact  that  long- 
time-average atmospherically  degraded  Images  are  extremely  difficult  to 
restore  In  practice.  In  principle.  If  enough  light  flux  were  utilized 
in  the  detection  of  the  Image,  significant  restoration  would  be  possibU 
However,  the  theoretical  results  above  imply  that,  for  the  conditlori; 
interest  here,  the  amount  of  light  flux  required  is  prohibitively  large. 

We  note  In  closing  that,  although  the  results  above  have  been 
derived  for  the  case  of  a po'  it-source  object,  they  can  be  shown  to  be  a 
close  approximation  for  the  case  of  a more  complicated  object  provided 


Oq/ o - I and  N is  replaced  by  n , the  average  number  of  photoevents 


per  correlation  area  of  the  object. 

5.  FUTURE  WORK 

The  formalism  described  above  is  now  ready  for  application  to  several 
problems  important  for  compensated  imaging.  We  are  now  in  the  process 
of  deriving  comparable  results  for  the  case  of  a partially  compensated 
imaging  system,  modeling  the  residual  wavefront  errors  as  a gaussian 
phase  screen.  In  addition,  we  intend  to  explore  the  restorabi I i ty  of 
short-exposure  images  using  the  model  developed  here. 
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APPENDIX 

Our  goal  in  this  appendix  is  to  show  that  the  fourth-order 
moment  of  Eq.(24)  i educes  to  the  result  indicated  in  Eq.(25).  Thus  we 
wish  to  find 


r LI  N N N N 

{idiN  “ I I I ):  t 

^ ■*  n=l  m=l  p=l  q=l  |_ 

* ”Y<''n'Vm  ♦ Vp-V,)]| 


exp|-j2ii  j^v^(x^"x_  + x_-x_) 


m 


P «i 


(A-l) 


The  N terms  in  this  summation  can  be  placed  in  15  different  classes 
as  fol lows: 


(1) 

n=m=p=q 

N terms 

(2) 

n*m,  p=q,  n^p 

N(N-I)  terms 

(3) 

n=m,  pi^qi^n 

N(N-l)(N-2)  terms 

(4) 

n=p,  m=q,  n>*m 

N(N-l)  terms 

(5) 

n=p,  mj<qj<n 

N(N-l)(N-2)  terms 

(6) 

n=q,  m=p,  nj*m 

N(N-I)  terms 

(7) 

n»q,  m>*p#n 

N(N-l)(N-2)  terms 

(8) 

n=m=p,  n>*q 

N(N-I)  terms 

(9) 

n«*m*q , ni*p 

N(N-I)  terms 

(10) 

n»p»q , nj*m 

N(N-l)  terms 

(11) 

p»q»m,  n/m 

N(N-I)  terms 

(12) 

ni^mi*pi<q 

N(N-l)(N-2)(N-3)  terms 

(13) 

p»q,  n>*mj<p 

N(N-l)(N-2)  terms 

(U) 

m»q,  ni*mi*p 

N(N-I)(N-Z)  terms 

(15) 

m*p,  n|(«n}(q 

N(N-l)(N-2)  terms 

I 

s 


We  regard  the  rate  X(x,y)  of  the  process  as  a known,  deterministic 

function,  and  average  over  the  2N'<-i  random  variables  (x^yj)  , 

(x-,y<,)  , ...  (x„,y,.)  , N.  Noting  that,  for  a Poisson  random  variable 
. N N 


N . 

E^[n(N-I  ....  (N-k+l)j  * [V)]*"  * 

the  contr ibut ions  of  the  15  sets  of  terms  are: 
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(A-2) 


■^1 

'■■'S 


(1) 

(2) 

"(X) 

L^fX). 

l\ 

• i 

. i. 

, ! 

s- 

(3) 

’"(X)’ 

^|A(Vj^,Vy)  1^ 

1 

1 

i 

i 

} 

W 

A(2vj^,2vy) 

i 

s 

(5) 

■"(X)] 

r A 

^A(2vj^,2vy)l  A (vjj,Vy) 

2 

(6) 

L"(x). 

i’_ 

1 

\ 

i 

/ 

(7) 

K)] 

^iA(Vj^.Vy)  1^ 

1 

(8) 

."(X). 

^lA(VjjfVy)  1^ 

} 

(9) 

'si 

^|A(Vj(,Vy) 

00) 

^lA(VjjtVy)  1 ^ 

j 

i 

(H) 

'si 

|^|A(Vjj»Vy) 

i 

02) 

.“(X)] 

^(A(Vjj,Vy)  1^ 

) 

5 

03) 

'"(X)] 

^|A(Vj5,Vy)  1^ 

i 

{ 

04) 

1 iv  r ^ * 

■^A  (2vjj»2vy)  1 A(vjj.Vy) 

2 

05) 

>(*)] 

^lA(WjjiVy) 

{ 

Here,  as  before,  the  definition 


n-j2it(v„x+v„y) 
X(x,y)  e dxdy 


A(vjj,Vy) 


* f 

X(x,y)dxdy 
j < 


is  used.  Noting  further  that 


A(vjj,Vy)  = N^^jA(vjj,Vy)  , 

and  combining  all  of  these  results,  we  obtain  Eq.(2S), 

+ |A(2Vjj,2Vy)|^ 

+ A(2Vjj,2Vy)  ^A*(vjj,Vy)j^  + A"(2vj^,2vy)  j^y\ 
+ |A(Vy,v„) . 
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I 3.  Single  image  estimate  of  spectral  density  illustrating  "half- 

I frequency"  noise,  (a)  sinusoidal  fringe  of  classical  intensity; 

(b)  single-image  estimate  of  the  spectral  density  of  that  intensity 
distribution. 
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4.  Block  diagram  for  least-mean'Square’'error  restoration. 
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(b) 


I 


5.  Normalized  spectral  densities  for  two  object  models: 

(a)  N independent  point  sources,  and  (b)  stationary  object  over 
a finite  region. 
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LONG  EXPOSURE 


6.  Maximum  restorable  frequency  AnZ/iT  vs.  N for  ■ 5 cm  and  10  cm, 

-C 

D - 152  cm  and  X - 5>‘10  ^ cm;  (a)  long  exposure,  no  tilt  removal, 

(b)  long  exposure,  tilt  removed,  far-field  atmospheric  propagation, 

(c)  long  exposure,  tilt  removed,  near-field  atmospheric  propagation. 
The  solid  horizontal  line  represents  the  telescope  cutoff  frequency 
0/X  , while  the  dashed  horizontal  lines  represent  r^/X  for  the  tv 
values  of  r^  used  here. 


